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Abstract 

Recent work by Kopp and MannhariP on novel electronic systems formed at oxide interfaces 
has shown interesting effects on the capacitances of these devices. We employ inhomogeneous 
dynamical mean-field theory to calculate the capacitance of multilayered nanostructures. These 
multilayered nanostructures are composed of semi-infinite metallic leads coupled via a strongly 
correlated dielectric barrier region. The barrier region can be adjusted from a metallic regime to a 
Mott insulator through adjusting the interaction strength. We examine the effects of varying the 
barrier width, temperature, potential difference, screening length, and chemical potential. We find 
that the interaction strength has a relatively strong effect on the capacitance, while the potential 
and temperature show weaker dependence. 



1 



I. INTRODUCTION 



As the capabilities of electronic components increase and their sizes decrease, new ideas 
are needed to continue the advancements in technology^. One of the fundamental electronic 
components is the capacitor, in its basic form it is a dielectric (k) layer separating two con- 
ducting layers, which builds and stores charge on each conducting plane when an external 
potential is applied. Driving the development of high-performance capacitors is the further 
miniaturization of various electronic devices including metal-oxide semiconductor field-effect 
transistors (MOSFETs). In order to continue the trend of development, high-K insulating 
materials are utilized in MOSFETs and capacitor s^"^. As the size of these devices decrease, 
quantum mechanical effects play a greater role and complicate the trend of using higher and 
higher k materials. Recent theory work by Kopp and MannharlP and subsequent experimen- 
tal work by Li et aP has introduced the idea of using ultrathin strongly correlated electronic 
materials to produce controllable small or large capacitances instead of the traditional high- 
K dielectric approach. This work builds on the growing number of possible applications of 
oxide interface^^^, as oxides forming two dimensional electron gases at the interface present 
a strong candidate for capacitance enhancements. The experimental work- found a greater 
than 40% increase in the gate capacitance when the mobile electrons were nearly depleted 
in a LaAlOs/SiTiOs interface. This increase is attributed to a negative compressibility of 
the interface electron system. Motivated by this work, we want to theoretically investigate 
the strong correlation effects on capacitance when the barier is a Mott insulator. 

We focus on constructing theoretical nanostructure devices consisting of ballistic metal 
leads on both sides of a strongly correlated electron dielectric layer. Inhomogeneous dy- 
namical mean-field theor}^^ (IDMFT) the theoretical framework for this work, allows for 
the self-consistent calculation of the properties of such devices. We use the Falicov-Kimball 
modeP^ to govern the interaction and use a Potthoff-Noltin^^' technique for solving the 
IDMFT. We work in the static limit with no current flow, where all calculations can be 
carried out in equilibriunP^. 

The capacitance is calculated for various parameters, showing the strongest dependence 
on the interaction strength and weaker dependence on the temperature and applied potential. 
Two methods for calculating the capacitance are discussed in the paper, one based on the 
center of charge approach of Lang and Kohri^ (where one measures the total charge on the 
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capacitor) and the other on the voltage profiles through the capacitor due to MeacP^ (where 
one measures the voltage difference between the plates). 

The general organization of the rest of this paper is as follows; in Section II, we detail 
the mathematical formalism and numerical issues associated with the calculation of the 
capacitance from the IDMFT approach. In Section III, we present numerical results of 
the capacitance dependence on various parameters including temperature, thickness, and 
dielectric screening length. We summarize the work, discussing the results and future ideas 
in Section IV. 

II. FORMALISM 

The general equation for the capacitance^, C, of two electrodes possessing charges of Q 
and —Qi separated by voltage V is 

C = ^. (1) 

In this work, we are concerned with parallel plate capacitors; an arrangement of two- 
dimensional layers stacked in a sandwich configuration. Classically, the capacitance of 
two-plate capacitor is defined as 

0="^, (2) 

where is the relative dielectric constant of the material separating the plates, eo (k = e^) is 
the dielectric constant of the vacuum, A is the area of the plates, and d is the thickness of the 
dielectric. These equations assume that the charges sit on an idealized surface plane of zero 
thickness. This assumption does not hold in reality and the electron density distribution 
must be taken into account. Kohn and Lan^^ showed that the effective position of the lead 
surface, be calculated from a center of charge approach, 

center center 
a=— oo a=— oo 

where is the position of plane a in the z-direction, and pa is the charge density distribution 
on plane a (with a = — oo being the left most plane and a = center being the center of the 
barrier region). The sum ranges over one half of the capacitor only and z^^ is equivalently 
defined with sum from the center of the barrier to last plane on the right. Equation ^ is 
therefore modified and the capacitance per unit area of a two plate capacitor becomes 

CcoC ^O^r 



A {z^ - 4) ■ 
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(4) 



For the rest of the paper we will refer to Eq. ^ as the center of charge (CoC) capacitance. 

In addition to the assumption of an ideal surface charge, as the size of these devices 
approaches an ultrasmall regime Eq. (|2| will also begin to break down. For example, a 
device where the effects of a single electron plays a dominate role, the calculations need to 
include quantum mechanics and a classical approach will not give a complete description. 
Quantum-mechanical effects will play an important role even before the single electron limit 
is reached. If the device is thin enough that the electric field is screened over a significant 
portion of the barrier region, there can be a noticeable reduction in the capacitance. This 
effect was first observed by MeacP^ for thin films and gives us a modified equation?^ for a 
thin film parallel plate capacitance per unit area; 

where Va is the applied potential, Vl and Vr are the potentials at the left and right interfaces, 
respectively, with {Vr — Vl) being the potential difference across the dielectric. Written in 
this form the modification to the geometric capacitance can be seen as the second term in 
the braces. We refer to Eq. (|5| as the voltage profile (VP) capacitance for the remainder of 
this paper. A schematic representation of the two methods is shown in Fig. [Tj 

The Cyp method can be thought of as fixing the distance between the plates to the 
physical difference of the plate geometry and calculating the effective potentials at that 
separation. While the CcoC method sums the total charge, Q, for the two halves and 
calculates an effective distance between charges. Hence, the effective charge in the VP 
method will not equal the total charge in the CoC method, while the distance between 
the plates in the CoC method will not equal the distance between the plates in the VP 
method unless the two results coincide. Each method has its relative strengths and relation 
to experiments discussed later in this paper. 

Extracting physical properties from inhomogeneous multilayered nanostructures is made 
possible by utilizing the algorithm originally employed by Potthoff and Noltin^^ and later 
adapted by Freericks^^. This so called quantum zipper algorithm reduces the complexity of 
the system by stacking two-dimensional translationally invariant planes, thereby building 
the inhomogeneity in the longitudinal third dimension {z) only. The z-coordinate remains in 
real space, while the x and y coordinates are Fourier transformed to wavevectors kx and ky, 
respectively, forming a mixed basis. Greek letters (a,/3,7, ...) are used to denote the planar 
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FIG. 1. Schematics of the two methods of calculation of the capacitance. Panel (a) shows the 
center of charge (CoC) method. The parameter Zq is calculated by using the charge distribution to 
the left of center and is shown as a dashed line. This line is shifted away from the original interface 
location shown as a dotted line. Similarly Zq is calculated by summing the charge distribution 
to the right of center. Figure (b) shows the voltage profile method (VP). Vl is the value of the 
potential at the left interface and Vr the potential at the right interface. Va is the applied voltage 
which equals two times AV. 
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index of the x — y planes stacked in the z-direction. We are left with a quasi-one-dimensional 
problem for each two-dimensional band energy, that can represented tridiagonally in real 
space and solved via the quantum zipper algorithm. The many-body equations are iterated 
to achieve a self-consistent solution. This Potthoff-Nolting approach is used to extract the 
electronic charge on each two-dimensional plane via the Green's function^. The charges 
are then used in a classical calculation to find the potentials on each plane and then the local 
electrochemical potential. The approach is iterated until it reaches a steady state charge 
distribution. 

We begin the calculations by defining the governing Hamiltonian for our system. The 
Hamiltonian involves a hopping term for the electrons and an interaction term for the sites 
within the barrier region. For the interaction in our numerical calculations, we employ the 
Falicov-Kimball modeP^ which involves an interaction between spinless conduction electrons 
and spinless localized electrons. When the conduction electron hops onto a site already occu- 
pied by the localized electron, both electrons feel a mutual repulsion. When this correlation 
strength is large enough in the Falicov-Kimball model it has a Mott-like metal-insulator 
transition. Additionally the Falicov-Kimball model has the advantage of simplifying the 
solution of the IDMFT equations for this system. It turns out that the stabilization of the 
voltage profile under the iterative solution of the IDMFT equations is difficult, and to date 
has only been able to be achieved in numerical solutions that are extremely accurate, like the 
Falicov-Kimball model. In the second quantization formalism, the spinless Falicov-Kimball 
Hamiltoniarp21 is, 



a ij^plane a idplane 



(6) 



a i£plane ol i£plane 

where the first two terms of the Hamiltonian describe intraplane and interplane electron 
hopping, respectively, where cjj^j and c^i are creation and annihilation operators on plane a 
and site i, respectively. The third term describes the charge reconstruction that occurs due 
to an externally applied potential, fia is the chemical potential on plane a, Va is the potential 
energy on plane a due to the Coulomb interaction of the electronic charge reconstruction, 
and AV^ is the input applied potential on plane a. The last term is the interaction term 
where Ua represents the interaction strength on plane a, and Wai is a classical variable that 



equals one if there is a localized particle at site i on plane a and zero if there is no localized 
particle at site i on plane a. 

We will use Green's functions to solve the problem. The equilibrium Green's function, in 
real space, and imaginary time is defined by 

Ga/3..M = -(TrC„.(r)4.(0)), (7) 

for imaginary time r, where Tr represents the time ordering operator. The notation (X) 
denotes the trace, Tr exp(— /3if)X divided by the partition function Z= Tr exp(— and 
the operators are expressed in the Heisenberg representation X(r) = exp(rif) X exp(— rif), 
all with respect to the Hamiltonian H. To properly express the Green's functions for the 
Matsubara frequencies we use a Fourier transformation 

G^p{iujn)=T dre'^-^G^p^r). (8) 
Jo 

where T = is the temperature. 

To build our model, we need to solve for the local Green's function on each plane, which we 
do by employing the quantum zipper algorithmpSl, based on the Potthoff-Nolting formalism?^. 
When solving for the local Green's functions, we use the fermionic Matsubara frequencies, 
iUn = i'n:T{2n + 1). We start with the unperturbed equilibrium equation of motion (EOM), 
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where e^n = — 2^^ [cos fc^ + cos ky], fc" = (fc^,, ky, 0) is defined as the transverse momentum, 
5q,/3 is the Kronecker delta function, and Haiiujn) is the self energy on plane a. Note that 
from this point on we will use the simplification that the hopping matrix elements are equal 
to t for nearest neighbors, ta+ia = ta-ia = taij = t and vanish otherwise. Since the EOM 
has a tridiagonal form with respect to the spatial component z{a, /3) it can be solved with 
the so-called renormalized perturbation expansiorP^l. We solve the equation directly, for the 
P = a case via 

G^JiUn, fc") = — 



(10) 



We create left and right recursion relations, 



La-nii(^n; fc") = + /i^ - K + AV^ - S«_„(iw„) - + -y-, (11) 



and 



^2 

Ra+n{i(^n] A;") = lUn + /ia - K + AV„ - Sa+„(icU„) - CfcH + — -ttt (12) 

respectively to solve for the other values oi a ^ (3. We start these relationships with the 
bulk values (n — )■ ±00), which give us 

L_oo(«a;„; fc") = ^ 

±^ V[2u;„ + ^ia-V^ + - S_oo(^c^n) - efcii]' - (13) 



and 



2 

± ^ ^[zWn + yUa - + AK. - Soo(^Wn) " efcll]' " ■ (14) 

The signs in the previous two equations are chosen to yield an imaginary part less than zero 
for iojn lying in the upper half plane, and vice versa for iun lying in the lower half plane. 
The self-energies vanish for the ballistic metal leads used here. 

To get our final expression for the Green's function, we substitute the left and right 



equations, Eqs. (11) and (12) respectively, into Eq. (10) to yield 



GaaiiOJn', k^^) = " • (15) 

The local Green's functions on each plane can then be found by summing the Green's 
functions over the transverse momenta 

with 



being the 2D density of states (DOS), K is the complete elliptical integral of the first kind 
and a is the lattice constant which we set to 1 for our calculations. After calculating the local 
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Green's functions on each plane, we use Dyson's equation to define the effective medium for 
each plane, 

GQ^{iUJn) = G~^{iUn) + T.a{iU}n)- (18) 

The local Green's function for the a-th plane then satisfies 

Wi 



Ga{tUJn) = -^Z^. ^ TTT- + 



(19) 



where wi is the average fiUing of the localized particles. Finally, we use the new local Green's 
functions and Dyson's equation again to find the self-energy. 



(20) 



This forms the basic algorithm for dynamical mean-field theory, which we now augment to 
determine the capacitance. 

To calculate the capacitance, we need the quantum-mechanically calculated electron num- 
ber density and the potential at plane a . Wc calculate the electronic charge on each plane 
by summing the Green's functions over all Matsubara frequencies on the imaginary axis, 
multiplied by the temperature. The electron number density at plane a satisfies 



Pa = 7; + T^Ga{iuJn)- 



(21) 



We can take advantage of the behavior of GaiioOn) at large n, which goes like l/iun, al- 
lowing us to regularize the Matsubara frequency summation by adding and subtracting 
^/[^^n + P — R^^a{ii^nmax)]- This givcs an exact summation of the tail of the Mat- 
subara sums and the electron number density becomes 

1 



Ga{iuJn) 



— - tanh 
2 



iun + Pa- ReEa(ia;n^^,) 
/3[Pa - ReEc,(ia;n^^J] 



(22) 



To find the Coulomb potential on each plane we begin with the magnitude of the local 
electric field created on plane a, 



lEl 



|e| \Po 



Jbulk I „ 
Pa " 



(23) 



2eoera 

where e is the charge of an electron, eo is the permittivity of free space, era is the relative 
permittivity of plane a, and p^"'*' is the bulk electron density of the material of which plane 



a is composed. Once the total field is known for each plane, we integrate them to find the 
electric potentials. Since the electric field's magnitude is constant, it is straightforward to 
compute the Coulomb potential. 



a 

where we define the parameter. 



E^J=a+i ^ [eschotil) + eschotil - 1)] , l3 > a 

0, 13 = a (24) 

E^""^ 2 i'^schotil) + eschotil + 1)] , l3 <a 



eschot{a) = , (25) 

which characterizes the decay of the surplus charge density away from the interface. The 
parameter p = (p„ — p''^''^)/N (with the total number of self-consistent planes used 
in our calculations) is used to improve the convergence of our equations and vanishes for 
the converged final fixed-point solution. It is worth noting that we are fixing eschot (equal 
in both the metal leads and barrier) and not recalculating the dielectric constant, meaning 
the many-body effects on the dielectric are already incorporated a priori in the calculation. 
The input eschot incorporates all contributions to the dielectric, including the bare dielectric, 
ion core, etc., therefore we can not directly compare our results to the geometric capacitance 
in Eq. ^ because we cannot isolate the different contributions to the dielectric to find the 
effective er needed in the formula for the geometric capacitance. 

We can now state the full algorithm used in our calculations. We start by inputting a 
value for the screening length (eschot), the applied potential {AV = Va/2), the chemical 
potential (/i), and the temperature (T). We begin the iterative calculations with a guess for 
the self-energy on each plane, usually zero or the solution to a previous calculation. Next, 
we use the left and right recursive equations to calculate the local Green's functions at each 
plane. These local Green's functions are then used to calculate the effective medium for each 
plane, which in turn is used to solve for the impurity Green's functions. The new impurity 
Green's functions are used to calculate the new self-energies which are used to feed the loop 
again. Additionally the impurity Green's functions are used to extract the planar filling. 
The planar filling is used within classical electrostatics to calculate the electric potential 
on each plane and in turn the contribution of the potential energy to the electrochemical 
potential on each plane. We average the new potentials, 

ynext iteration ^ ^^yold ^ _ ay)V^ew ^26) 
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with a large damping factor ay- The parameter ay is usually at least 0.99 which is needed 
to slow the updating and prevent converging to a nonphysical solution. We iterate through 
these steps until the calculations converge. Due to the large damping factor, it typically 
takes between 1,000 and 10,000 iterations to reach convergence. In order to achieve proper 
convergence the errors are kept to less than one part in 10^, this allows for reproducibility 
in the algorithm ensuring the planar charge densities are accurately calculated from one 
iteration to the next. Keeping the error below a tolerable level has proven difficult to 
achieve in models other than the Falicov-Kimball model (such as the Hubbard model). 

After the calculations reached a self consistent solution we can use our two methods to 
extract the capacitance. The CoC capacitance is calculated by summing over the extracted 
planar filling [Eq.(|4]) for both the right and left halves of the layers]. The VP capacitance 
is calculated by substituting the calculated electric potential at the left and right interface 
layers in for Vl and Vr, respectively [Eq. (Is])]. 



III. RESULTS 

All of our numerical results will be calculated at half-filling (/i = 0, (c|cj) = 1/2, and 
Wi = {wi) = 1/2). We build our model with 30 self-consistent metal planes in the leads 
each terminating in the bulk surrounding the dielectric layers in the center. We vary the 
thickness of the dielectric region from 4 to 20 planes. The calculations are carried out on a 
simple cubic lattice allowing only nearest neighbor hopping (both interplane and intraplane 
hopping, t, are equal). This reduces the number of parameters in the calculations allowing 
focus on the properties of interest. 

There are many parameters that can be varied to investigate their effects on the ca- 
pacitance. For each set of parameters in our calculations, we can extract the quantum- 



mechanically calculated electron number density for each plane from Eq. (22), as seen in 



Fig. |2] (a), which plots the difference between the electron number density and the bulk 



electron number density through the device for various Cschot- From Eq. (24) we can also 
plot the potentials on each plane through the device, which is shown in Fig. [2] (b), again for 
various values of eschot- Note how the most rapid change in the potentials occur near the 
interface. With the charge and the potentials known, we can calculate the capacitance for 
each device. Once the potentials have been calculated through the devices the capacitance 
is calculated from Eq. ([s]). Figure [2] (a) and Figure |2] (b) show that as the screening length 
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increases the charge deviation curves become, as expected, sharper in nature. 

Once we have calculated the capacitance per unit area we can compare the center of 
charge capacitance as defined in Eq. ^ with the voltage profile capacitance calculated 
from Eq. ([s]). Both methods of calculating the capacitance are plotted in Fig. [s] which 
shows that as we vary the screening length, both formulas for the capacitance follow the 
expected C/A oc (l/eschot) behavior. 

We investigate thermal effects on the capacitance by varying the temperature. Figure 
|4] shows that the capacitance/area falls off slightly (~ 4%) as the temperature is increased 
over the range of 0.1 to 0.8. In our calculation, if we take a reasonable energy scale for 
our system, such as a noninteracting bandwidth of 3 eV, then t = 0.25eV {U = 6). This 
corresponds to a temperature range from room temperature (T = 0.1) to 2500K (T = 0.8). 

Fig. [5] shows that the capacitance per unit area has a stronger dependence on the Falicov- 
Kimball interaction strength (U). The CoC capacitance per unit area grows faster than the 
VP capacitance per unit area. The CoC capacitance per unit area increases by approximately 
30% across the range of interaction strengths. There is a crossover around U = 7.5 where 
the center of CoC becomes larger than the VP capacitance. The increase in capacitance that 
results from increasing the interaction strength in the barrier is due to the higher interaction 
strength reducing the spatial extent of the dipole layer that is formed on the inside of the 
dielectric region. 

We plot the charge deviation and potential profiles for various applied potentials {V) in 
Fig. |6]as well as the capacitance/area for various applied potentials in Fig. [Tj The reduction 
in capacitance due to increasing the applied potential can be seen in Fig. |6] as the planes on 
either side of the dielectric layer begin to saturate and excess charge is forced further away 
from the interface leading to a reduction in the capacitance. 

IV. CONCLUSIONS 

We presented in this paper a self-consistent method for using IDMFT to calculate the 
capacitance of multilayered nanostructures. We also discussed the various capabilities and 
challenges with this many-body approach. We showed how using the quantum zipper al- 
gorithm based on the work of Potthoff and Nolting, we can calculate the electron number 
densities on each plane. The electron number densities were then used to calculate the 
electric potential on each plane from classical electrostatics. We presented two methods to 
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FIG. 2. (Color online) [10 dielectric planes, U = 6, wi = 0.5, Ej = 0, K = |2|, and T = 0.25] 
Charge deviations from the bulk charge (a) and the electric potential (b) plotted as a function of 
plane numbers for various eschot as indicated by the legends. 
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FIG. 3. (Color online) [10 dielectric planes, U = 6,wi = 0.5, Ef = 0, = |2|, and T = 0.25] The 
capacitance per unit area plotted as a function of eschot- 
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FIG. 4. (Color online) [10 dielectric planes, eschot = 1, U = 6, wi = 0.5, Ef = 0, and Va = |1|] 
The capacitance per unit area plotted as a function of temperature. 
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FIG. 5. (Color online) [10 dielectric planes, eschot = 0.5, wi = 0.5, Ef = 0, Va = and 
T = 0.25] Capacitance per unit area for various interaction strengths U, with a crossover seen at 
approximately U = 7.5 



calculate the capacitance, one based on a center of charge approach and the other account- 
ing for the screening of the charges. The capacitance is calculated for various parameters, 
finding the strongest dependence on the interaction strength (U), 30% over the calculated 
range. We find a weaker dependence on temperature (5-10%) and applied potential. 

We use a semiclassical approach to calculate the potential, which fixes the relation be- 
tween the applied potential and the electronic charge density at a given plane. By making 
the dielectric a fixed parameter we cannot capture any effects on the capacitance that are 
not already described in Eq. ^ or Eq. (|5]). Fixing the dielectric in the system forces 
the parameter to describe the total dielectric not just the dielectric values coming from the 
polarizability of the ion cores. Although this limits our ability to isolate many-body effects 
in Eq. ^ or Eq. (|5| from the geometric capacitance in Eq. (|2]), we are still able to examine 
the behavior of the model as we vary other parameters. 

We calculated the capacitance via two methods, through the generated voltage profiles, 
Cyp, as well as the center of charge approach, CcoC- Calculating the capacitance in these 

15 




30 40 

Plane Number 



CO 

c 
0) 

o 

CL 

"to 
.g 

'l_ 
-I— « 
O 
0) 

LU 






-aV=0.5 


— ■ 


-aV=2 




-AV=3.5 


— •- 


-aV=6 


—k- 


-aV=10 


— ▼- 


-AV=14 



Plane Number 



FIG. 6. (Color online) [10 dielectric planes, eschot = 1, f/ = 6, = 0.5, Ef = 0, and T = 0.25] 
Charge deviation from the bulk (a) and the electrical potential (b) profile for various applied 
potentials. 
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FIG. 7. (Color online) [10 dielectric planes, eschot = 1, f/ = 6, = 0.5, Ef = 0, and T = 0.25] 
Capacitance/ Area for various applied potentials, Va- 



methods allows for comparison to different experimental setups. For example, if the capac- 
itance of an experimental set up is determined by integrating the total charge required to 
discharge the capacitor, then comparison to the CcoC method would be more appropriate 
because the total charge is pa and we are not measuring the potential at the capacitor 
plates but between the left and right leads. On the other hand, if the experimental set 
up had the ability to probe the potential at the edge of the metallic plates, then the Cyp 
method is more appropriate for comparison, since we do not measure the total charge. 

The ability to calculate the potential and charge profiles that incorporate many-body 
effects allow us to investigate any non-linear effects. The calculations for the capacitance 
presented in this paper were carried out in the slow limit with no current floW^, reproducing 
the inverse dependence on thickness and linear dependence on eschot- In future research 
moving away from half-filling would allow for the model to hopefully capture more non- 
linear behavior as the Falicov-Kimball model enters into a phase-separated state. Other 
potentially interesting and non-linear behavior can be investigated by using the Hubbard 
model rather than the Falicov-Kimball model, which is currently inhibited by the achievable 
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accuracy of the required numerical calculations. 

The experiments see the greatest enhancements to the capacitance near the highly de- 
pleted limit, therefore by moving the Falicov- Kimball model away from half filling and 
forcing the model to enter a phase-separated state, we expect to see enhancements in our 
model. The many body effects should become more pronounced as a phase-separated state 
is entered. 
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